Purpose: compare normalization methods for a particular dataset

source(file.path("scripts", "CZI_functions.R"))
if (!("Rtsne" %in% installed.packages())) {
   install.packages("Rtsne")
}
if (!("NMI" %in% installed.packages())) {
   install.packages("NMI")
}
if (!("caret" %in% installed.packages())) {
   install.packages("caret")
}
# Get the file names of all the normalized files
normalized.files <- dir("darmanis_data/normalized_darmanis")

# Read in each of the normalization files
normalized.data <- lapply(normalized.files, function(file) {
                              readr::read_tsv(file.path("darmanis_data/normalized_darmanis",
                                                        file))
                          })
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   gene = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   Genes = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   Genes = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   Genes = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   Genes = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   Genes = col_character()
## )
## See spec(...) for full column specifications.
# Keep all the gene lists
genes <- lapply(normalized.data, function(x) x[,1])

# Keep the names with it.
names(normalized.data) <- gsub("\\.tab|darmanis_", "", normalized.files)
# Read in the data
meta <- readr::read_tsv(file.path("darmanis_data", "GSE84465_meta.tsv"))
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   channel_count = col_double(),
##   taxid_ch1 = col_double(),
##   contact_phone = col_double(),
##   `contact_zip/postal_code` = col_double(),
##   data_row_count = col_double(),
##   `plate id:ch1` = col_double(),
##   `tsne cluster:ch1` = col_double()
## )
## See spec(...) for full column specifications.
# Get sample names from columns
samples <- colnames(normalized.data[[1]])[-1]

# Keep metadata only for the samples we have
meta <- meta[match(samples, meta$geo_accession), ]

# Extra cell types info as it's own vector
cell.types <- as.factor(meta$`cell type:ch1`)
plate.batch <- as.factor(gsub("plate id: ", "", meta$`characteristics_ch1.1`))
# Run tsne on each dataset and extract
tsne <- lapply(normalized.data, function(dat) {
                        tsne.res <- Rtsne::Rtsne(t(dat[,-1]),
                                                 check_duplicates = FALSE)
                        tsne.res <- tsne.res$Y
                        names(tsne.res) <- samples
                        
                        return(tsne.res)
                        })
# Make a function for plotting tsne's and using 
tsne.plot <- function(dat, var, name = "name"){
  # This function is to plot tsne data and label it by a variable
  # Args:
  #  dat: a data.frame with two columns of data
  #  var: vector contains metadata labels
  #  name: a character string for the plot title and for to save as png
  # Returns:
  #   png and plot in Rmd of the provided data with labels of the metadata provided
  # convert celltype into numbers
  colz <- colors(distinct = TRUE)[runif(length(levels(var)), min = 1,
                                        max = length(colors(distinct=TRUE)))]
  plot(dat,pch = 21, bg = colz[var], main = name);
  legend(x = "bottomleft", legend = levels(var), fill = colz, cex = 0.6)
}

# Plot em with cell type and plate batch labels
lapply(tsne, function(dataset) {
      # Get normalization method
      set.name <- names(tsne)[parent.frame()$i[]]
      
      # Make plots for cell type and plate batch 
      cell.type.plot <- tsne.plot(dataset, cell.types, name = set.name)
      plate.batch.plot <- tsne.plot(dataset, plate.batch, name = set.name)
    
      # Save the plots to pngs
      png(paste0("tsne_", set.name, "cell_type.png"))
      cell.type.plot
      dev.off()
      
      png(paste0("tsne_", set.name, "plate_batch.png"))
      plate.batch.plot
      dev.off()
      
      # Print out the plots in the Rmd
      cell.type.plot
      plate.batch.plot
})

## $counts
## $counts$rect
## $counts$rect$w
## [1] 8.141351
## 
## $counts$rect$h
## [1] 94.51097
## 
## $counts$rect$left
## [1] -24.61422
## 
## $counts$rect$top
## [1] 71.69175
## 
## 
## $counts$text
## $counts$text$x
##  [1] -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946
##  [8] -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946
## [15] -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946
## [22] -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946
## [29] -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946
## [36] -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946
## [43] -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946 -22.20946
## [50] -22.20946
## 
## $counts$text$y
##  [1]  69.8385972  67.9854409  66.1322846  64.2791283  62.4259719
##  [6]  60.5728156  58.7196593  56.8665030  55.0133467  53.1601904
## [11]  51.3070341  49.4538778  47.6007215  45.7475652  43.8944089
## [16]  42.0412526  40.1880962  38.3349399  36.4817836  34.6286273
## [21]  32.7754710  30.9223147  29.0691584  27.2160021  25.3628458
## [26]  23.5096895  21.6565332  19.8033769  17.9502205  16.0970642
## [31]  14.2439079  12.3907516  10.5375953   8.6844390   6.8312827
## [36]   4.9781264   3.1249701   1.2718138  -0.5813425  -2.4344988
## [41]  -4.2876552  -6.1408115  -7.9939678  -9.8471241 -11.7002804
## [46] -13.5534367 -15.4065930 -17.2597493 -19.1129056 -20.9660619
## 
## 
## 
## $DESeq2
## $DESeq2$rect
## $DESeq2$rect$w
## [1] 6.615268
## 
## $DESeq2$rect$h
## [1] 112.4473
## 
## $DESeq2$rect$left
## [1] -25.8007
## 
## $DESeq2$rect$top
## [1] 84.78113
## 
## 
## $DESeq2$text
## $DESeq2$text$x
##  [1] -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671
##  [8] -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671
## [15] -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671
## [22] -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671
## [29] -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671
## [36] -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671
## [43] -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671 -23.84671
## [50] -23.84671
## 
## $DESeq2$text$y
##  [1]  82.5762771  80.3714288  78.1665804  75.9617321  73.7568838
##  [6]  71.5520355  69.3471872  67.1423389  64.9374906  62.7326423
## [11]  60.5277940  58.3229457  56.1180974  53.9132491  51.7084008
## [16]  49.5035525  47.2987042  45.0938559  42.8890076  40.6841593
## [21]  38.4793110  36.2744627  34.0696144  31.8647661  29.6599178
## [26]  27.4550695  25.2502212  23.0453729  20.8405246  18.6356763
## [31]  16.4308280  14.2259797  12.0211314   9.8162831   7.6114348
## [36]   5.4065865   3.2017382   0.9968899  -1.2079584  -3.4128067
## [41]  -5.6176550  -7.8225033 -10.0273516 -12.2321999 -14.4370482
## [46] -16.6418965 -18.8467448 -21.0515931 -23.2564414 -25.4612897
## 
## 
## 
## $log2
## $log2$rect
## $log2$rect$w
## [1] 8.306228
## 
## $log2$rect$h
## [1] 90.09418
## 
## $log2$rect$left
## [1] -26.15157
## 
## $log2$rect$top
## [1] 64.74442
## 
## 
## $log2$text
## $log2$text$x
##  [1] -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812
##  [8] -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812
## [15] -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812
## [22] -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812
## [29] -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812
## [36] -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812
## [43] -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812 -23.69812
## [50] -23.69812
## 
## $log2$text$y
##  [1]  62.9778700  61.2113175  59.4447650  57.6782125  55.9116600
##  [6]  54.1451075  52.3785550  50.6120025  48.8454500  47.0788975
## [11]  45.3123450  43.5457925  41.7792400  40.0126875  38.2461350
## [16]  36.4795825  34.7130300  32.9464775  31.1799250  29.4133725
## [21]  27.6468200  25.8802675  24.1137150  22.3471625  20.5806100
## [26]  18.8140575  17.0475049  15.2809524  13.5143999  11.7478474
## [31]   9.9812949   8.2147424   6.4481899   4.6816374   2.9150849
## [36]   1.1485324  -0.6180201  -2.3845726  -4.1511251  -5.9176776
## [41]  -7.6842301  -9.4507826 -11.2173351 -12.9838876 -14.7504401
## [46] -16.5169926 -18.2835451 -20.0500976 -21.8166501 -23.5832026
## 
## 
## 
## $scVLM
## $scVLM$rect
## $scVLM$rect$w
## [1] 8.859139
## 
## $scVLM$rect$h
## [1] 84.06793
## 
## $scVLM$rect$left
## [1] -31.52045
## 
## $scVLM$rect$top
## [1] 58.87337
## 
## 
## $scVLM$text
## $scVLM$text$x
##  [1] -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368
##  [8] -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368
## [15] -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368
## [22] -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368
## [29] -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368
## [36] -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368
## [43] -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368 -28.90368
## [50] -28.90368
## 
## $scVLM$text$y
##  [1]  57.2249841  55.5765932  53.9282024  52.2798115  50.6314206
##  [6]  48.9830297  47.3346388  45.6862480  44.0378571  42.3894662
## [11]  40.7410753  39.0926845  37.4442936  35.7959027  34.1475118
## [16]  32.4991209  30.8507301  29.2023392  27.5539483  25.9055574
## [21]  24.2571666  22.6087757  20.9603848  19.3119939  17.6636030
## [26]  16.0152122  14.3668213  12.7184304  11.0700395   9.4216487
## [31]   7.7732578   6.1248669   4.4764760   2.8280851   1.1796943
## [36]  -0.4686966  -2.1170875  -3.7654784  -5.4138693  -7.0622601
## [41]  -8.7106510 -10.3590419 -12.0074328 -13.6558236 -15.3042145
## [46] -16.9526054 -18.6009963 -20.2493872 -21.8977780 -23.5461689
## 
## 
## 
## $TMM
## $TMM$rect
## $TMM$rect$w
## [1] 9.576485
## 
## $TMM$rect$h
## [1] 91.43549
## 
## $TMM$rect$left
## [1] -33.33695
## 
## $TMM$rect$top
## [1] 65.31656
## 
## 
## $TMM$text
## $TMM$text$x
##  [1] -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829
##  [8] -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829
## [15] -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829
## [22] -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829
## [29] -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829
## [36] -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829
## [43] -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829 -30.50829
## [50] -30.50829
## 
## $TMM$text$y
##  [1]  63.5237050  61.7308521  59.9379993  58.1451465  56.3522937
##  [6]  54.5594409  52.7665881  50.9737352  49.1808824  47.3880296
## [11]  45.5951768  43.8023240  42.0094711  40.2166183  38.4237655
## [16]  36.6309127  34.8380599  33.0452071  31.2523542  29.4595014
## [21]  27.6666486  25.8737958  24.0809430  22.2880902  20.4952373
## [26]  18.7023845  16.9095317  15.1166789  13.3238261  11.5309733
## [31]   9.7381204   7.9452676   6.1524148   4.3595620   2.5667092
## [36]   0.7738563  -1.0189965  -2.8118493  -4.6047021  -6.3975549
## [41]  -8.1904077  -9.9832606 -11.7761134 -13.5689662 -15.3618190
## [46] -17.1546718 -18.9475246 -20.7403775 -22.5332303 -24.3260831
## 
## 
## 
## $voom
## $voom$rect
## $voom$rect$w
## [1] 8.410643
## 
## $voom$rect$h
## [1] 83.74553
## 
## $voom$rect$left
## [1] -27.94895
## 
## $voom$rect$top
## [1] 58.74276
## 
## 
## $voom$text
## $voom$text$x
##  [1] -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465
##  [8] -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465
## [15] -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465
## [22] -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465
## [29] -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465
## [36] -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465
## [43] -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465 -25.46465
## [50] -25.46465
## 
## $voom$text$y
##  [1]  57.1006897  55.4586204  53.8165511  52.1744818  50.5324125
##  [6]  48.8903432  47.2482739  45.6062046  43.9641353  42.3220659
## [11]  40.6799966  39.0379273  37.3958580  35.7537887  34.1117194
## [16]  32.4696501  30.8275808  29.1855115  27.5434422  25.9013728
## [21]  24.2593035  22.6172342  20.9751649  19.3330956  17.6910263
## [26]  16.0489570  14.4068877  12.7648184  11.1227491   9.4806797
## [31]   7.8386104   6.1965411   4.5544718   2.9124025   1.2703332
## [36]  -0.3717361  -2.0138054  -3.6558747  -5.2979440  -6.9400134
## [41]  -8.5820827 -10.2241520 -11.8662213 -13.5082906 -15.1503599
## [46] -16.7924292 -18.4344985 -20.0765678 -21.7186371 -23.3607065
kmeans_eval <- function(feature, metadata = metadata, iter = 10, seed =1234) {
  # This function is used to perform iterative k-means clustering based on projected 
  # features for single cell data and then evaluate the performance according to 
  # Nomalized mutual information (NMI) and adjusted rand index (ARI)
  #
  # Args:
  #  feature: a data.frame contains projected features (n dimensional space, n = 2, 3 ...), the columns are
  #    projected features in n dimensional space for a sample (cell), rows are samples
  #  celltype: vector contains cell type informantion
  #  iter: number of interation for k-means clustering
  #  seed: seed for k-means clustering
  # Returns:
  #   NMI and ARI results
  # convert celltype into numbers
  metadata_num <- as.numeric(factor(metadata))
  sample_id <- seq(1:nrow(feature))
  
  # iterative k-means
  nmi_score_all <- c()
  ari_score_all <- c()
  all_cluster <- list()
  
  for(i in 1:iter){
    # set k equal to the number of celltypes in the dataset
    k <- length(unique(metadata))
    
    # perform k means clustering
    km <- kmeans(feature, k)

    # true clusters
    orignal_data <- data.frame(sample_id, metadata_num)
    
    # predicted clusters
    cl_data <- data.frame(sample_id, km$cluster)
    
    # calculate NMI and ARI score
    nmi_score <- NMI::NMI(orignal_data, cl_data)$value
    ari_score <- mclust::adjustedRandIndex(km$cluster, metadata_num)
    
    nmi_score_all <- c(nmi_score_all, nmi_score)
    ari_score_all <- c(ari_score_all, ari_score)
  }

  # Compile all results into a data.frame
  results <- data.frame(ari = ari_score_all, nmi = nmi_score_all)
  return(results)
}
knn_eval <- function(feature, metadata = metadata, k = 10){
  # This function performs knn based evaluation 
  # Args:
  #  feature: a data.frame contains projected features (n dimensional space, n = 2, 3 ...), the columns are
  #    projected features in n dimensional space for a sample (cell), rows are samples
  #  celltype: vector contains cell type information
  #  k: cross validation fold
  # Returns:
  #  list of accuracy scores for each iteration of KNN
  # Make the data into a data.frame:
  feature <- data.frame("tsne" = feature, "metadata" = metadata)
  
  # Split observations into groups
  cv <- cvTools::cvFolds(nrow(feature), K = k, R = 1)
  
  # Create empty objects to store the performance information for each iteration
  perf.eval <- list()
  confusion.matrix <- 0
  
  # Go through this iteration that k times
  for (i in 1:k) {
    # Isolate samples for training the model
    train <- feature[cv$subsets[-which(cv$which == i)], ]
    
    # Isolate samples for testing the model
    test <- feature[cv$subsets[which(cv$which == i)], ]
    
    # Perform KNN model fitting
    knn.fit <- caret::train(metadata~. , data = train, method = "knn",
                     trControl = caret::trainControl(method = "cv", number = 3),
                     preProcess = c("center", "scale"),
                     tuneLength = 10)
    
    # Evaluate the model
    knn.pred <- predict(knn.fit, newdata = subset(test, select = -c(metadata)))
    perf.eval[[i]] <- round(cal_performance(knn.pred, test$metadata, 3), 2)
    
    # Make the results into a matrix
    matrix <- as.matrix(table(test$metadata, knn.pred, deparse.level = 0))
    confusion.matrix <- matrix + confusion.matrix
  }
  
  # Get mean performance of cross validation
  perf.eval <- dplyr::bind_rows(perf.eval)
  accuracy <- perf.eval$accuracy

  return(data.frame("knn" = accuracy))
}
# Get knn and kmeans results for all tsne's of all datasets
cell.type.results <- lapply(tsne, function(dataset) {
                        knn.results <- knn_eval(dataset, metadata = cell.types)
                        kmeans.results <- kmeans_eval(dataset, metadata = cell.types)
                            
                        data.frame(knn.results, kmeans.results)
                            })
## Loading required package: lattice
## Loading required package: ggplot2
# Get knn and kmeans results for all tsne's of all datasets
batch.results <- lapply(tsne, function(dataset) {
                        # knn.results <- knn_eval(dataset, metadata = plate.batch)
                        kmeans.results <- kmeans_eval(dataset, metadata = plate.batch)
                            
                        # data.frame(knn.results, kmeans.results)
                        data.frame(kmeans.results)

                            })
plot.results <- function(results.list, name = "results") {
  # This function makes a boxplot of the cluster statistic results 
  # Args:
  #  results.list: a list of dataframes which each column contains a different 
  #                statistic for 10 iterations
  #  name: name to use for the png to be saved and the plot title 
  # Returns: boxplots of the normalization method cluster statistics. Prints the
  #          plots and save the plot as a png
  # Get meta info 
  meta <- names(unlist(results.list))

  # Transform list into a dataframe:
  ggplot.df <- data.frame("method" = stringr::word(meta, 1, sep = "\\."),
                          "test" = gsub("[0-9]*$", "", 
                                               stringr::word(meta, 2, sep = "\\.")),
                          "iter" = rep(1:nrow(results.list[[1]]),
                                       length(results.list)*ncol(results.list[[1]])),
                          "values" = unlist(results.list))
  
  # Make the plot
  plot <- ggplot(data = ggplot.df, aes(x = method, y = values, fill = test)) +
          geom_boxplot(position = position_dodge()) + 
          xlab("Normalization method") +
          ggtitle(name) +
          facet_wrap(~test)

  # Save plot to png
  ggsave(paste0(name, "cluster_results.png"), width = 10) 
  
  # Print plot
  plot
}
# Plot the cell type results 
plot.results(cell.type.results, name = "Cell_type")
## Saving 10 x 5 in image

# Plot the plate batch results
plot.results(batch.results, name = "Plate_batch")
## Saving 10 x 5 in image

sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin17.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libLAPACK.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] caret_6.0-81    ggplot2_3.1.0   lattice_0.20-35
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_0.2.5   purrr_0.2.5        reshape2_1.4.3    
##  [4] splines_3.5.1      colorspace_1.3-2   generics_0.0.2    
##  [7] stats4_3.5.1       htmltools_0.3.6    yaml_2.2.0        
## [10] prodlim_2018.04.18 survival_2.42-6    rlang_0.3.0.1     
## [13] e1071_1.7-0        ModelMetrics_1.2.2 pillar_1.3.0      
## [16] NMI_2.0            withr_2.1.2        glue_1.3.0        
## [19] bindrcpp_0.2.2     foreach_1.4.4      plyr_1.8.4        
## [22] bindr_0.1.1        lava_1.6.4         robustbase_0.93-3 
## [25] stringr_1.3.1      timeDate_3043.102  munsell_0.5.0     
## [28] gtable_0.2.0       recipes_0.1.4      codetools_0.2-15  
## [31] evaluate_0.12      labeling_0.3       knitr_1.20        
## [34] class_7.3-14       DEoptimR_1.0-8     Rcpp_0.12.19      
## [37] readr_1.3.0        scales_1.0.0       backports_1.1.2   
## [40] ipred_0.9-8        hms_0.4.2          digest_0.6.18     
## [43] stringi_1.2.4      Rtsne_0.15         dplyr_0.7.7       
## [46] grid_3.5.1         rprojroot_1.3-2    tools_3.5.1       
## [49] magrittr_1.5       lazyeval_0.2.1     tibble_1.4.2      
## [52] crayon_1.3.4       cvTools_0.3.2      pkgconfig_2.0.2   
## [55] MASS_7.3-51        Matrix_1.2-14      data.table_1.11.8 
## [58] lubridate_1.7.4    gower_0.1.2        assertthat_0.2.0  
## [61] rmarkdown_1.10     iterators_1.0.10   mclust_5.4.1      
## [64] R6_2.3.0           rpart_4.1-13       nnet_7.3-12       
## [67] nlme_3.1-137       compiler_3.5.1